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THE  GENERALIZED  INVERSE  PROBLEM  OP  ORBIT  COMPUTATION  « 


V.  O.  SZEBBHELY  . 

Space  Sciences  Laboratory,  ^Qtneral  Electric  Company,  Philadelphia,  Pa.,  USAJ 

Abstract;  The  basic  problem  of  celestial  mechanics  —  finding  the  motion  for  given 
initial  conditions  in  a  given  force  field  —  is  inverted  by  asking  for  force  fields 
which  result  in  either  a  specifio  orbit  or  which  allow  to  obtain  a  closed  form  solution 
of  the  equations  of  motion.  The  first  case  lesuls  to  a  generalized  concept  of  guidance 
and  thnut  programs,  while  the  second  opens  the  way  to  finding  new  closed  form 
solutions  for  sligiitly  modihod  n-body  (n  >  2)  gravitational  fields. 

The  problem  discussed  in  this  paper  is  a  generalization  in  purpoe)  and  in  method 
of  the  classical  problem  of  orbit  determination  in  celestial  mcohanics.  It  allows 
the  determination  of  classes  of  functions  which  either  furnish  optimization  oppor¬ 
tunities  of  the  guidance  problem  (first  case)  or  furnish  overall  field  modifications 
which  allow  closed  form  solutions  (second  case)  applicable  to  perturbation  cal¬ 
culations  as  reference  orbits  with  predetermined  approximations. 

Elstablishing  new  reference  orbits  for  perturbations  and  formulating  a  concise 
and  systematic  method  of  guidance  analysis  are  the  two  main  results  presented 
in  the  paper. 

PeawMe:  Oniicunaorcn  MeroA,  c  noMOutbio  Koreporo  onpexMiniOTCir  MOXiKtm- 
leaitHH  cun,  AeflcTnyioutMX  Ha  KOCMHqecKHe  puKeru  tskhm  oCpaaoM,  ito  hx 
opCiiTu  H  TpacKTopiiM  MoryT  Curb  oiiitcaiiiJ  m  npeACTaDAciiu  c  noMotubio 
pcuiciiiiii  li  saMKiivToll  ({lopMe  c  AaitituM  npHflAHHceiiHCM.  OcHODiiaii  aoAaqa 
HoOeciiofl  McxaiiMKH  -  naxoHCAeitite  ADH}icei]HJi  xan  aaiihux  HaioJibHUX 
ycnooHll  n  AaiiuoM  ciiaooom  hoao  -  oOpamoeTCH  M  cuoamtch  k  saAaqe  iiaxoHC- 
ACHHfl  ciiAouoro  noAH,  KOTopoe  xacT  n  peayAbTare  hah  onpcAe-ieiiiiyio 
opOHTy,  HAH  noanoAncT  noAymiTb  peuieime  n  aaMKiiyroll  ijiopMe  ypaniieitHil 
ABHHcenHH.  B  nepnoM  cAynae  AaiiHoe  ciiAonoe  hoao  MOAHiJiHiutpyercJi  oxoxb 
TpaeKTopiiH,  HTopaH  aaAamt  TpeCyer  MOAH<|)HKauHH  noAiioro  hoah.  IlepBufl 
CAynail  ueAer  k  oCoCmeiiHoft  KOimenuHii  nporpaMM  ynpaoACHHH  h  thth, 
TOPAa  KaK  HTopoll  oTKpuDaeT  nyxii  k  iiaxoHtACHHio  iiodux  peuieiiHfl  b  aaii- 
Kiiyrofl  (|iopMfl  aah  caaGo  MOAH(|)HUHpbDniiiihix  rpaBHTauHOHHiJX  noAell  n 
TCAa  (n  >  2).  IvAaccimecKaii  aaAaMo  onpeAOAeiiHfl  opCiiru  b  iieOecHoll  iiexa- 
ilHKC  COCTOHT  B  HaXO>KACHHH  aACMCHTOn  OpGliTU  H  (HAH)  liaHaAbllUX  yCAOBHfl 
iia  ocHoue  iiaOAioAeHHii  n  chaobom  hoac,  T.e.  npcACTanAJieT  coCofl  aaAaxy 
oOpariiyK)  aaAanc  pacnera  rpaeKTopiiH  aah  Aaiinoro  chaobopo  hoah  h  xniiiiba 
HaMOAbiibix  ycAOBiiil.  Qua  uo)kot  Curb  iioAH<|)HitHpouaHa  c  ueAbio  iiaxoMC^ 
ACHHH  OnpeABACHIIMX  KOHCTUHT  HAH  XUpUKTCpHCTHK  CHAOBOPO  HOAH,  SaKOH 
THPOTCiiHH  KOTopopo  npcAnoAapaeTCfl  aiiAaiiiibiu  (iianpiiMcp,  noAo  SeuAH 
HAH  oTKpiJTKe  IlonTyiia).  PnccMaTpHsaeMAH  n  crarbe  oCpaTiioH  saxasa 
npcAcrauAHOT  coGuli  oGoGiueiiHe  ho  hoah  h  Meroxy,  tuk  kak  ohu  nosBOJiHeT 
onpcAOAHTb  K-mccu  (|)yHKmiil,  KOTopiJO  oGocncMHnaioT  ahGo  onTHMaAbiiue  Boa- 
MOHtiiocTH  iipoGACMU  ynpaBACiiHH  (nepnijfl  CAyHafl),  ahGo  oGimte  moah(|)H- 
KaitHH  HOAH,  noanoAHiouHiu  noAyMHTb  pcuieiiHJi  b  sauKHyToA  i^iopMe  (nropofl 
cAyMaii),  npHMeiiHMiJo  k  pacMoraM  c  BoaMyuteHHiiMH  hcxoahux  opGht  c 
aapaiiee  onpoACAeiiHMMH  iiphCaiijkahhmmh.  ycTaiioBAeiiHe  iiooux  hcxoahmx 
opGiiT  iipit  iinJiiiqiii!  BoaMyiueiiiin  ii  (fiopMyAHponaiiHe  Tomioro  n  cHcreMa- 
TimecKOPO  MCTOAa  aiioAHaa  HABeAeiiHH  -  abb  ochobiimx  pesyAbTara  HacTOHiqefl 
puCoTU. 

*  This  publication  stems  in  part  from  work  sponsored  by  tbs  United  States  Air 
Force,  contract  AF-49(638)-814,  Air  Force  Office  of  Scientific  Research  of  the  ARDC. 
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1.  Introduction 

Inverse  techniques  are  widely  used  in  several  branches  of  theoretical 
physics,  applied  mathematics,  continuum  mechanics,  celestial  mechanics,  etc. 
In  some  of  these  fields  the  inverse  approach  has  lead  to  the  solution  of 
important  problems,  in  others  it  opened  up  entirely  new  fields  of  research 
and  in  a  few  cases  the  potentialities  have  not  yet  been  fully  evaluated.  The 
more  flexible  the  concept  of  inverse  technique  is,  the  more  usefulness  it 
offers,  therefore,  in  the  following,  a  general  view  will  be  presented  and 
applied  to  some  basic  problems  of  applied  celestial  mechanics. 

A  phenomenon  is  often  described  by  a  set  of  diffeiential  equations  aii^ 
initial  and  boundary  conditions.  A  “solution”  is  often  pre  ented  by  a  set  of 
independent  variables  as  functions  of  space,  time  or  other  dependent  variables. 
These  functions  u'e  considered  solutions  if  they  satisfy  the  above  mentioned 
equations  and  conditions.  The  idea  of  the  inverse  approach  stated  in  its  most 
general  form  is  to  start  out  with  the  “solution”  and  find  the  equations 
describing  the  phenomenon.  Inasmuch  as  this  general  statement  is  of  very 
limited  use  for  solving  actual  problems,  certain  restrictions  are  applied.  In 
hydrodynamics  for  instance,  the  continuity  equation  for  incompressible 
potential  flow  is  the  Laplace  equation  which  the  potential  function  (<p)  must 
satii^fy.  Therefore,  functions  which  satisfy  the  /i(p==0  equation  will  describe 
certain  flow  fields.  The  real  and  imaginary  parts  of  functions  of  a  complex 
variable  satisfy  the  Laplace  equation,  therefore,  flow  fields  m’ght  be  generated 
inversely  by  finding  the  boundary  conditions  which  are  satisfied  by  arbitrary 
chosen  functions.  To  satisfy  part  of  the  problem  by  giving  solutions  of  the 
differential  equations  and  then  finding  the  boundary  conditions  is  a  well 
known  approach  in  continuum  mechanics. 

On  the  other  hand,  to  find  the  equation  which  will  bo  satisfied  by  a  “solu¬ 
tion”  which  in  turn  satisfies  certain  initial  conditions  is  a  frequently  en¬ 
countered  problem  in  celestial  mechanics.  Even  the  above,  already  re  tricted 
—  as  compared  to  the  original  —  statement  is  too  general,  since  often  only 
certain  constants  in  the  differential  equations  aie  found  by  impo.sing  a  solution. 
Frequently  the  general  form  of  the  equations  are  given  and  the  “solution” 
will  be  selected  so  that  it  will  further  spcH-ify  the  equations.  The  classical 
example  is,  of  com^e,  when  observational  data  represent  the  “solution”  and 
the  equations  describing  the  phenomenon  are  sought.  Some  of  the  physical 
laws  governing  the  phenomenon  are  known  and  some  ate  unknown,  some 
numerical  constants  attached  to  the  problem  are  given,  some  aie  to  be 
determined.  Classical  examples  are  the  discovery  of  Neptune  from  its  pertur¬ 
bative  effect  on  the  orbit  of  Uranus,  establishment  of  the  coefficients  of  the 
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'  higher  order  harmonics  in  the  expansion  of  the  Earth's  gravitational  poten¬ 

tial,  etc. 

Another  example  of  somewhat  different  type  is  the  solution  of  the  earth 
satellite  problem  by  Vinti  [1]  and  Gaifinkel  [2].  Slight  modifications  of  the 
potential  function  describing  the  Earth’s  gravitational  field  either  by  adjusting 
numerical  coefficients  or  by  changing  the  functions  —  might  result  in  satisfy¬ 
ing  the  new  equations  by  a  “solution”  which  would  not  satisfy  the  original 
equations. 

Several  questioas  remain  open  in  connection  with  the  above  mentioned 
examples.  Differential  equations  describe  physical  phenomena  only  approx¬ 
imately  and  it  can  hap])en  that  a  certain  set  of  data  will  “fit”  a  solution  of 
the  modified  equations  better  than  it  satisfies  the  original  equations.  If  the 

I  numerical  constants  occurring  in  an  equation  are  modified  and  this  way  a 

solution  is  obtained,  the  deviation  must  be  evaluated  between  the  solutions 
to  the  original  and  to  the  modified  equations.  If  not  only  numerical  constants 
but  the  functions  involved  are  modified,  this  evaluation  might  be  very 
difficult. 

The  avenues  of  ap[)lications  of  the  inverse  approach  in  celestial  mechanics 
have  boon  varied  and  proven  to  be  extremely  useful.  The  purpose  of  this 
paper  is  to  present  a  general  view  and  t<o  point  out  its  applications. 

2.  Analytical  considerations 
2.1.  Elimination  of  time 

Consider  dymamical  systems  of  n  degrees  of  freedom  with  generalized 
coordinates  q\,  qz . qu  ...,  qn,  and  write  the  equations  of  motion  in  the  form: 

=  (1) 

Considering  eveiy  qt  depending  on  qi,  we  have 

doi 

and 

or 

(2) 

where  the  prime  denotes  derivatives  with  respect  to  qi  and  the  dot  denotes 
derivatives  with  respect  to  the  time  (t).  From  (2): 

For  t»2,  eqs.  (1)  and  (3)  give 

ftiqk,  7»)+7i'(?i)* 


(3) 
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or  using  eq.  (2): 


hiqk,  qk,  qi)=q2fi{qk.  qk,  gi)+g2'(?i)*. 

Solving  this  equation  for  qi,  or  for  (gi)*,  we  write 

g/,  qt) 

and  computing  the  time  derivative,  we  obtain 

Eqs.  (1),  (2),  and  (6)  transform  (6)  into 


k-i  ^k^^^  k.ti^qk  d/  '^6^2*  d/ 


S21 

(4) 

(6) 

(6) 


Oi  /  /  •v  m 

1  55;  <»  +  1  »•  +  5^  »  <’) 

where  91'  se  1,  provided  that  qi  ^  0. 


For  t  — 3,  eq.  (3)  gives 

hiqk,  qk,  qt") = qtfxiqk,  qk,  qk’) + qt’Hiqt,  qt,  gz’).  (8) 

The  form  of  the  equation  is  the  same  for  t>3(l)n: 

IMk,  qk,  qt') = qi’hiqk,  qk,  qt') + qt'Higk,  qk,  qt')-  (0) 


Eq.  (7)  is  a  third  order  differential  equation  and  the  (n-2)  eqs.  (9)  are 
of  the  second  order.  The  system  is  of  (2n—  1)  -th  order  as  expected  since  it 
was  obtained  by  elimination  of  the  time  from  the  original  (2n-th  order) 
set.  It  is  significant  that  eq.  (7)  contains  only  one  third  derivative.  In  the 
two  d^rees  of  freedom  case  1,  2  and  only  eq.  (7)  is  obtained  by  elimination 
of  the  time.  This  equation,  however,  will  describe  completely  the  geometry 
of  the  problem. 

Let  q\-x,  qt  =  y,  and  i  =  fiix,y,x,y)  ) 

y=h{x,y,i,y)S 

oorrosponding  to  eqs.  (1).  The  simple  relations  of  y  =  y'x  with  y’=dy/dx 
and  y=y(x),  correspond  to  eqs.  (2).  Eq.  (4)  becomes: 

ft(x,  y,  y',  i) = y7i(*.  y.  y'.  + y'(*)*  ( i  * ) 

which,  when  solved  for  x  for  (x)*  gives: 

(x)*  =  ff(x,y,y',y').  (1*) 


Final  elimination  of  the  time  results  in  the  equation  of  the  orbit,  corres- 


322 

ponding  to  eq.  (7); 


V.  O.  SZEBEHELY 


[ra] 


The  above  equatioas  of  motion  (10)  represent  a  slightly  more  general  case 
than  the  two  dimensional  restricted  three  body  problem  in  a  rotating  rectan¬ 
gular  coordinate  system.  The  differential  equation  connecting  y  and  x  is  of 
the  third  order.  If  uSi)  is  made  of  the  Jacobi  integral,  eq.  (13)  can  be  reduced 
to  the  second  order,  as  will  be  shown  later.  The  three  dimensional  case  of 
the  restricted  three  body  problem  corresponds  to 

x  =  h{x,y,y)  I 

y  =  f2{x,y,x)  /  (14) 

2  =  fsix,  y)  ) 


and  it  can  be  shoum  that  the  6-th  order  system  is  reduced  to  the  5-th  order, 
consisting  of  a  third  order  and  of  a  second  order  differential  equation.  The 
first  equation  contains 

X,  y,  z,  y',  z',  y%  z\  y" 

and  the  second 

X,  y,  2,  y\  z\  y',  z\ 

The  following  three  remarks  are  offered  at  this  point: 

(i)  By  eliminating  the  time  two  purposes  were  served:  (a)  the  order  of 
the  system  of  differential  equations  has  been  reduced  by  one,  and  (b) 
since  time  effects  do  not  enter,  the  problem  is  changed  into  a  geo¬ 
metrical  om. 

(ii)  If  the  equations  of  motion  contain  the  time  explicitly,  eq.  (6)  will 
contain  an  additional  term  dH'bt,  in  which  case  elimination  of  time  is 
not  achieved. 

(iii)  Since  the  elimination  of  time  process  implies  that  the  qi  =  qi(t)  function 
be  inverted  to  t  =  t(qi),  and  .substituted  in  the  ?j  =  ?j(0  relation,  giving 
qi  =  qi[t(qi)]  the  (/i  =  0  |)oint  requests  siHXjial  attention.  Double  valued 
functions  ocf’urring  in  eq.  (6)  warrant  the  same  comment. 

2.2  Modified  equations 

Consider  the  right  side  of  eqs.  (1)  as  the  sum  of  two  functions: 

fi(qk,  qn)  =  Ftiqt,  qk)  +  fiiqk) 

it,  i  =  l(l)n 


(16) 
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i.e  the  modifying  fimctions  {<pi)  depend  only  on  the  generalized  position 
coordinates.  The  effect  of  the  set  will  now  be  investigated  on  the  geometry 
of  the  orbit. 

Eq.  (4)  will  assume  the  form  of 

Fiiqit,  qk,  qi)+<pt(qk)=qz[Fi{qk,  qt,  ?i)+9i(?*)]+?a'(?i)*.  (16) 

From  this  equation  we  obtain  91,  assuming  that  F\{qk,  qt,  qi)  and 
Ft{qk,qk,qi)  ftre  given  and  the  q>i(qk)  and  q>2{qk)  functions  are  not. 

qi*==H{qk,  qk  ,  qz",  fi,  q>2)-  (17) 


The  first  equation  of  the  orbit  corresponding  to  eq.  (7)  is  obtained  again 
by  differentiation: 

nr,  «  ,  •  bH  ,  6H  „  bH  •  b<pi  ,  ] 

2Fi  +  2<pi  2  r—  qk  +  ^  -pT—,  qk  +  jr— ;  ?2  +  r —  J  t—  qk  j 

fc-i  o^t  k-z^qk  bqz  k-i 

For  given  Fi,  Ft  and  q)\,  q>2  functions,  eq.  (18)  is  one  of  the  orbit  equations. 
For  given  Fi,  Ft  and  qk^^qkiqi)  functions,  eq.  (18)  is  a  partial  differential 
equation  for  <pi  and  <pt.  For  given  Fi  and  Ft,  eq.  (18)  is  an  orbit  equation 
and  can  be  used  to  select  and  (pt  to  facilitate  obtaining  the  solution. 
The  other  orbit  equations  corresponding  to  eq.  (9)  become  for  »=s3(l)n: 

Fi  +  (pi  —  qt{Fi  +  fpi)i^qt''H  (19) 

i.e.  every  orbit  equation  introduces  one  new  pt  function. 

The  previously  mentioned  set  of  eqs.  (10)  in  a  modified  form  becomes: 

x  =  Fi(x,y,x,y)  +  fi(x,y)}^ 
y-=r:{x,  y,  x,  y)  +  <ft(^,  y).i 
Eq.  (16)  for  this  case  is 

Ft  +  <pt=y'{Fi  +  pi)  +  y''(x)^  (21) 

from  which  x*  can  be  computed,  giving 

x^=H{x,y.y',y'',pi,<pt).  (22) 

The  equation  of  the  orbit,  corresponding  to  eq.  (18)  is 

„„  „  bH  ,bH  ,bH  „bU  \ 

2f.+2,i-'5j+y5j+»  5^+y  J 


\  dz  ^  by  /  bq>i  \bx  ^  by  J  bq>2 


(23) 
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This  is  MU  ordinary  third  order  ditferential  equation  for  t/(z),  i.e.  for  the 
function  representing  the  orbit,  if  Fi,  Ft,  (pi  and  (ft  are  given.  It  also  can  be 
considered  as  a  partial  differential  equation  for  <p\(x,  y)  and  ipt(x,  y),  i.e.  for 
the  modifying  functions,  if  Fi,  Ft  and  y=y{x)  are  given.  Finally,  eq.  (23) 
can  also  serve  as  a  guide  in  selecting  the  fi,  and  <pt  functions  for  given  Fi 
and  Ft  so  that  the  solution  can  be  represented  by  functions  of  known  geo¬ 
metric  properties. 

As  summary  of  this  chapter,  the  following  remarks  are  offered: 

(i)  For  an  n  degree  of  freedom  system  there  are  n  modifying  functions: 
q>i...<Pn  which  depend  on  the  n  generalized  position  coordinates, 
qi  ...  qn.  After  elimination  of  the  time  the  2n-th  order  system  is  reduced 
to  a  system  of  (2n—  l)-th  order,  represented  by  n— 1  equations.  The 
first  equation  is  of  the  3-rd  order.  It  includes  (pi,  tpt  and  their  partial 
derivatives,  the  latter  ones  linearly.  The  remaining  (n  — 2)  equations 
each  include  971,  (pt  and  only  one  of  the  ft  functions  {k=Z{l)n). 

(ii)  If  the  Ft  functions  and  the  qt(qi)  functions  are  given  the  n  modifying 
functions  ft  will  have  to  be  determined  from  n  —  1  equations,  one  of 
these  being  a  partial  differential  equation,  i.e.  the  ft  functions  are  not 
uniquely  determined. 

(iii)  For  given  Ft  functions  (properly)  selected  ft  functions  will  completely 
determine  the  solution. 

3.  Applications 

3.1.  General  considerations 

It  will  be  shown  that  the  above  conclusion  marked  (ii)  allows  the  formu¬ 
lation  of  generalized  guidance  equations  on  one  hand  (section  3.1.1)  and 
offers  a  new  presentation  of  the  Encke  method  on  the  other  hand  (section 
3.1.2).  Conclusion  (iii)  points  the  way  to  the  establishment  of  new  approximate 
solutions  (section  3.1.3). 

3.1.1.  Oeneralized  guidance  approach 

The  combination  of  eqs.  (1)  and  (16)  describes  a  dynamical  system  where 
the  “forcing  functions”  are  split  into  two  parts: 

qt  =  F liqt,  qt)  +  ftiqt)-  (24) 

The  first  part  (Ft)  is  now  considered  as  the  given  force  field,  the  second 
part  (ft)  as  the  guidance  force  along  an  orbit.  To  determine  the  n  guidance 
force  components  (ft),  eqs.  (18)  and  (19)  are  used.  Since  there  are  only  (n- 1) 
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equations,  no  unique  solution  is  available,  in  fact,  the  guidance  forces 
fi.  <PA  •••  fn  will  be  expressed  as  functions  of  <p\  and  <pt.  It  is  noted  that  fp\ 
and  (ft  are  not  independent  but  related  by  a  partial  differential  equation. 
From  the  point  of  view  of  gi<:idance  this  is  one  of  the  advantages  of  the  present 
approach  since  a  great  variety  of  guidance  functions  result  from  the  analysis. 
In  actual  cases,  the  ftrst  guidance  function  might  be  determined  as  a  function 
of  the  position  coordinates  and  of  the  second  guidance  function  (902)-  Even 
if  ipt  is  fixed,  the  first  guidance  function  is  not  determined  uniquely  since 
it  appears  as  a  solution  of  a  partial  differentia)  equation.  This  equation  has 
the  general  form: 


(26) 


where  L  is  a  linear  function  of  tlie  partial  derivativos. 

The  determination  of  the  set  of  guidance  functions  therefore  is  redticeil 
to  first  finding  the  general  solution  of  eq.  (25),  for  instance,  in  the  form  o 
—  Using  then  oqs.  (19)  the  remaining  guidance  functions  are 

determined  as 


<Pi=‘fiisPi^  qi)^ 
« =  3(  1  )n.  S 


(26) 


Optimizing  conditions  imposed  on  the  problem  might  require  finding  the 
minimum  of  the  total  force  vector,  {Wi)  of  work  {W2)  or  of  power  {Wa),  i.e. : 


=  +  ...  or  'P2  =  i^<pkdqk  or  ^  ftdqk.  (27) 

3.1.2.  Oeneralized  Encke  method 

The  Encke  perturbation  method  computes  the  difference  between  the 
actual  and  a  fictitious  precomputed  (“nominal”)  trajectory.  The  coordinates 
of  this  nominal  trajectory  can  be  iepre.sented  by  closed  form  solutions  of  a 
differential  equation.  The  method  is  based  on  the  premi.se  that  if  the  differen¬ 
tial  equation  lor  the  nominal  trajectory  is  a  “slight”  modification  of  the 
ifferential  equation  associated  with  the  original  jiroblem,  then  the  deviation 
between  the  actual  and  the  nominal  trajectories  is  small. 

The  Encke  method  [3]  is  summarize!!  here  for  comparison  using  the 
notations  of  this  pajier  and  projicrly  .selected  units: 


let  qt=  -  ^  ^ 

r' 


(28) 


where  qi  is  the  t-th  [po.sition  coordinate,  qt^  and  Pt  is  the  t-th  com¬ 

ponent  of  the  perturbation.  These  equ^tioas  represent  the  actual  trajectory 
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while  the  reference  trajectory  satisfies  the  following  differential  equation 


ro® 


(29) 


where  p(  is  the  i-th  coordinate  of  the  two  body  reference  orbit  and 

rf?=  2 

4-1 

Introducing  for  the  difference  between  the  actual  and  the  reference  trajec¬ 
tories  ;ri  one  obtains  the  differential  equation,  which  the  Encke  method 
integrates : 


7*4  = 


(30) 


714  =  9*  —  />4.  ^ 

For  actual  computational  purposes  eq.  (30)  is  written  as 

..a  *  .  o 

ro*  fo® 


(31) 


where  the  computation  of 


a=  1  — 


[1(7*4 +  P4)*/'' 

is  facilitated  by  tables. 

If  the  reference  orbit  is  “close”  to  the  actual,  large  integration  steps  can 
be  taken,  therefore,  the  efficiency  of  the  Encke  method  depends  on  the 
proper  choice  of  the  reference  orbit.  This  is  not  difficult  when  planetary 
motions  are  studied,  i.e.  when  frequent  rectifications  (selection  of  new 
reference  orbits)  are  not  necessary.  For  purposes  of  lunar  trajectories,  or 
in  general,  for  three  body  orbit  calculations  the  original  Encke  scheme  has 
limited  significance  especially  in  the  regions  where  the  de>'iations  from  two 
body  orbits  are  largo. 

It  seems  to  be  a  logical  extension  of  the  Encke  method  to  investigate 
reference  orbits  which  are  not  based  on  two  body  calculations  but  which 
approximate  the  trajectories  associated  with  the  three  body  fields.  Let 

qi^tMk,qk)  (32) 

be  the  i-th  differential  equation  describing  the  actual  motion  of  a  body  in  a 
given  (/<)  force  field  and  let  the  reference  orbit  be  described  by 

pi=-fi{pk,  Pk)  +  fi(Pk) 

where  p<  is  the  i-th  coordinate  of  the  reference  orbit. 


(33) 
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The  partial  differential  equations  for  (tpi)  which  are  associated  with  eq.  (33) 
are  obtained  from  eqs.  (18)  and  (19)  by  writing  for  qt  and  /<  for  Ff  -. 

2/i  +  2^1  =  XI  (pfc.  Pk\  Vk,  Pa",  Vu  f2,  ^  (34) 

and 

/<  +  9’<  =  P<'(/i  +  9’i)  +  P<'f/(p*,  P*',  Pa*,  fi,  9>2),  for  i  =  3(l)n.  (36) 

A  set  of  simple  p<(pi)  functions  can  be  constnicted  which  represent 
approximately  the  solution.  Substituting  these  functions,  eqs.  (34)  and  (35) 
will  give  several  sets  of  <pi  functions.  These  sets  are  subjected  to  a  minimiza¬ 
tion  process  and  a  single  set  of  (pt  functions  is  obtained.  This  set  of  q>t  with 
the  associated  p<(pi)  functions  complete  the  constiuction  of  the  leference  orbit 

Introducing  now  the  deviations  of  the  actual  orbit  from  the  reference 
orbit  by 

7ii  —  qt  —  Pi 

and  subtracting  eq.  (33)  from  (32)  one  obtains; 

iii  =  fi(gk,  qk)-fi{pk,  Pk)-fi(Pk)  (36) 

where  jidO  is  to  be  determined,  =  P*(0  and  pt(l)  is  known. 

This  concludes  the  generalization  of  the  Encke  method  excepting  the 
important  fact  that  Encke  writes  the  final  eq.  (31)  in  a  form  sjiecifically 
suited  for  numerical  work.  The  corresjionding  step  in  this  generalized  treat¬ 
ment  can  be  accomplished  only  if  the  sjiecific  /(  functions  describing  the  field 
of  an  actual  problem  are  given. 

3.1.3.  Approximate,  solutions 

The  problem  of  approximate  solutions  is  twofold:  firstly,  a  method  of 
constructing  such  solutions  is  to  be  established  and  secondly,  the  accuracy 
is  to  be  estimated.  It  will  be  assumed  that  the  solution  of  the  original  system 
of  differential  equations  is  not  available,  so  the  estimation  of  the  error  of 
the  approximate  solution  must  be  more  sophisticated  than  comparing  two 
functions  in  a  given  domaiii.  Error  estimates  are  intimately  associated  with 
the  functions  (9?d  previously  called  guidance  or  modifying  functions  as  was 
shown  before. 

In  section  3.1.1  the  guidance  problem  was  discussed  and  it  w'as  shown 
how  the  fi  functions  are  determined  along  a  prescribed  trajectory,  i.e.  how 
the  given  force  field  is  modified  along  a  given  orbit.  The  present  chapter 
deals  with  the  modification  of  the  entire  force  field,  so  that  approximate 
orbits  subject  to  arbitrary  initial  conditions  can  be  represented  in  “closed” 
form. 


? 
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The  Encke  method  is  a  special  |)orturbation  technique  and  the  aim  of  the 
previous  section  (3.1.2)  was  to  generalize  its  reference  orbit  aspects.  The 
purpose  of  the  present  chapter  is  to  show  how  the  inverse  trajectory  concept 
can  be  applio<l  to  obtain  approximate  general  solutions  and  how  to  estimate 
the  accuracy  of  these  solutions. 

Eqs.  (18)  and  (19)  with  —  0  repre.sent  the  orbit  in  a  given  Ft  force 
field  and  it  is  assumed  that  the  qi  =  qi(qi)  functions  which  satisfy  these 
equations  with  given  initial  conditions  are  not  available.  In  order  to  find 
an  approximate  solution  to  the  9:^  =  0  actual  problem,  one  finds  the  (pt  ^  0 
sets  for  which  the  general  solution  of  the  orbit  differential  equations  can  be 
given.  In  other  words,  the  exact  general  solution  for  the  orbit  in  a  modified 
field  Fi  +  is  determined  first.  The  approximate  general  solution  of  the 
actual  problem  is  represented,  therefore,  by  the  exact  general  solution  of  the 
approximate  problem. 

Since  the  difference  betw'een  the  solutions  of  the  actual  and  modified 
differential  eqiiations  are  in  general  proportional  to  the  modifying  functions 
a|)plied,  the  motion  in  the  original  force  field  will  differ  from  the  motion  in 
the  modified  field  by  an  amount  which  is  proportional  to  the  field  distortion. 
The  (ft  functions  rejiresent  the  distortion  of  the  field  and  the  smaller  these 
distortions  are,  the  better  the  approximation  becomes. 

7’he  description  cf  the  analysis  involved  is  simple.  Referring  to  eqs.  (18) 
and  (19)  one  sele<'ts  an  arbitrary'  set  of  9:1(9*)  functions  for  which  the  91(91) 
solutions  can  lie  repre.sented  by  “simple”  functions.  This  does  not  require 
the  preci.se  specification  of  the  qn  set,  only  the  general  functional  forms  of 
its  members.  The  next  step  is  to  minimize  the  values  of  the  members  of  the 
qi  field  distortion  set  in  the  domain  of  interest.  The  result  will  be  a  set  o 
qi  functions  which  will  allow  a  simple  representation  of  the  general  solution 
(91)  and  at  the  same  time  will  minimize  the  deviation  from  the  actual  solution. 

Three  remarks  are  of  some  importance: 

(i)  The  error  estimate  for  the  approximate  general  solution  is  associated 
with  a  domain  in  which  the  solution  is  applicable. 

(ii)  The  error  estimate  in  general  will  be  very  crude  unless  the  special  form 
of  the  differential  equations  is  taken  into  account. 

(iii)  The  field  distortion  functions  might  introduce  large  deviations  if  the 
domain  is  “large”.  Approximate  solutions  consisting  of  several  functions 
which  are  matched  at  the  boundaries  of  their  respective  domain  might 
have  to  be  considered. 
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3.2.  The  planar  re.stricted  three  body  problem 
3.2.1.  Reduction  to  the  second  order 


The  differential  equation.s  of  motion  of  the  title  jiroblein  in  a  rotating 
Cartesian  rectangular  coordinate  system  with  origin  at  the  ma.ss  center  of  the 
two  principal  bodies  (see  fig.  1)  can  bo  found  in  any  standard  reference  [4] 
in  the  following  form; 


£  =  n^x  +  2ny  — 


v(x-b) 

pS 


y  =  n^y-'2nx-  ^ 


(37) 


where  n  is  the  angidar  velocity  of  the  system,  r  is  the  distance  between  mi 
and  m,  q  is  the  distance  between  ni^  and  m,  =  v  =  k^m2 


a 


lm2 

ftii  T  m2 


and 


6  = 


Irtii 

mi  +  m2 


(  38 


Fig.  1.  C(X)nlinato  gyBtom  for  the  restricted  3-bo<ly  problem. 


The  Jacobi  integral  is 

(i)2  4  (y}2  -  n^x2  +  y2)-^  -  ^  =  2C.  (39) 

The  third  order  differential  equation  of  the  orbit,  containing  x,  y,  y',  y" 
and  y"  is  obtained  by  using  eq.  (7).  The  order  might  be  reduced  to  two  if 
use  is  made  of  the  Jacobi  integ  al.  These  steps  w'ill  bo  executed  below,  but 
first  a  transformation  to  polar  coordinates  will  bo  made,  by  the  relations; 

X  t  a  =  r  cos  0  and  y  =  r  sin  0. 


(40) 
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After  substitution,  the  equations  of  motion  (37)  become: 

o  /I  v{r  —  lcoH0)  \ 

r  =  r(<)  4  n)2  —  an2  cos  6— ^ - ; - •'  j 

r*  ^3  J 

) 


0=  -2  - (()  +  n)  + 


an*  sin  0  vl  sin  0 


r  r  rp-* 

ami  the  Jacobi  integral  assumes  the  form: 


(r)*  4-  (r6)2  _  jt2;-2  4_  2ran*  cos  0  —  —  —  —  =  2C 

r  0 


where 


It  is  noted  that  using 


p*  =  r*  —  2rl  cos  0  +  P. 


V  =  —  -  --  4  ratp  cos  0 
r  Q 


eqs.  (41)  and  (42)  can  be  written  as 

f  =  r(tf +  n)2-  ^ 
or 


rd 


and 


/ 

(r)*  +  (r<))2-n2r*  +  2K  =  2C. 


(41) 


(42) 

(43) 

(44) 


(46) 


(46) 


It  is  noted  that  the  (r,  0)  polar  coordinate  system  is  originated  at  mi  and 
not  at  the  origin  of  the  (x,  y)  coordinate  system  and  therefore  direct  comparison 
of  eqs.  (46)  and  (46)  with  reference  [4]  might  be  misleading.  The  third  order 
differential  equation,  corresponding  to  eq.  (7)  is  obtained  if  the  following 
substitutions  are  made; 


qi  =  0, 


qz^r,  /i=- 


2r0  +  n) 
r 


r*  be  ’ 


fz  =  r{0  +  n)2  - 


dr 


In  order  to  arrive  at  a  second  order  differential  equation,  use  will  be  made 
also  of  the  Jacobi  integral.  Writing  f  =  r'6  where  r'  =  (dr)j(dd),  (d)*  can  be 
expre.s.sed  using  the  Jacobi  integral  (46): 

2(C-F)  +  n*r*  ^  g*(r.  d) 

r*4(r')=  r*  +  (r')3  ^  ^ 

where  the  g{r,  0)  function  was  introduced  for  simplicity’s  sake  and  its  physical 
meaning  is  the  magnitude  of  the  velocity  of  the  mass  m  in  the  rotating 
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coordinate  8y8tem.  It  is  noted  that  ^  has  not  been  expressed  from  an  equation 
corresponding  to  (4)  in  the  general  treatment  since  the  Jacobi  integral  could 
be  used  instead. 

Corresponding  to  eq.  (3)  we  have 

f  =  r'04-r'{6)2  (48) 

or  according  to  (46) 

r(^  +  n)2  -  =  -  ^  [2r(<)  4  n)  ^  p  ^]  +  r'{6)K  (49) 

This  equation  corresponds  to  eq.  (4).  If  the  Jacobi  integral  would  not  exist, 
one  would  solve  eq.  (49)  for  (d)*  and  would  proceed  as  in  the  general  case. 
Using  the  expression  for  (6)*  as  given  by  eq.  (47).  observing  that 


dV 

dr 

'^  2^ 

bV 

60 

2  60 

9ir,0) 

T  —  1 

lr*  +  (r')2’ 

and  substituting  in  eq.  (49)  the  result  is: 

J  (r-  -  2  -  r)  +  [r>  +  (r’)"]  ^  [r^  +  (60) 


This  second  order  diflferential  equation  for  the  r  =  r(6)  function  represents 
the  orbit  associated  with  the  planar  restricted  three  body  problem. 

If  the  reciprocal  transformation  for  the  radius  vector  is  introduced, 
eq.  (60)  assumes  the  following  form: 

j7(m' +  «)  =  [^^^^  +  ^  l]"  (51) 

where  u=ljr,  the  subscripts  represent  partial  derivatives,  and 

fl2  fi  m2 

o*  =  2C  4  —  —  2 - cos  0  +  2m  [m  4-  v(  1  —  2lu  cos  0  +  I^m*)"’''*].  (62) 

u 


Eq.  (61)  is  the  result  of  the  derivation.  Its  general  solution  would  represent 
the  solution  of  the  restricted  three  body  problem. 
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Two  remarks  related  to  eq.  (51)  might  be  in  order; 

(i)  Due  to  the  choice  of  the  coordinate  system  and  because  of  the  structure 
of  the  differential  equation,  the  latter  reduces  easily  to  the  corresponding 
two  body  problem  by  writing 

r  =  0  and  n  -■  0, 

For  this  case,  which  corresponds  to  the  motion  of  tn  in  the  field  of  mi 
in  a  fixed  coordinate  system,  eq.  (52)  becomes: 

g^  =  2ftu4  2C  and  so  gu  =  filg-  (53) 

Eq.  (51)  gives: 

u’ +  u=[(uy^u^]filg^  (54) 

or 

u  +  =  (55) 

where  h  is  the  constant  of  integration  of  the  momentum,  i.e. 

A  =  dr*  =  6lu*. 


It  is  noted  that  eq.  (47),  using  the  (u,  0)  variables  can  be  written  as 


«))■- 


and  therefore 


g^u* 

u»4  (u')» 


(56) 


(u')*  +  u* 


which  relation  was  used  to  obtain  eq.  (55). 

The  fact  that  eq.  (54)  is  identical  with  the  well  known  form  of  the  two 
body  equation,  i.e. 

M  +  u'  =  constant 


can  also  be  shown  without  resorting  to  the  momentum  integral  as  follows: 

Eq.  (54)  can  be  written  itr 

/i[(uy  +  «*]  =  2(C  +  /4U)(U 4  u') 
which  by  differentiation  becomes 

/i[(u')*  ^  «*]'  =  2/iu'(u  4-  u')  4  2(C  4-  /iu)f u  4  u')'.  (57) 

Comparing  (57)  with  the  identity: 

[(tt')*+tt*]'  =  2tt'(u4-u')  (58) 

we  have  u  +  u' >=  constant,  provided  that  the  velocity  is  not  zero. 
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(ii)  The  Lagrangian  solutions  foUow  rather  elegantly  from  eq.  (61)  if 
u  =  uo  ^  constant  solutions  are  searched  for.  Since  the  g{u,  0)  function 
is  identical  with  the  velocity  magnitude  in  the  rotating  system  and 
since  the  Lagrangian  solution.^  are  stationary  in  this  system, 


g  =  0. 


This  reduces  the  differential  eq.  (61)  to 


/iWo®  H - -.{l—uol  cos  Oo)  =  n2(  1  —  nuo  cos  Oo). 

QfT 


(59) 

(60) 


Eliminating  a  and  by 

„._!L  and  (61) 

fi  +  v  F 

and  introducing  the  a  =  /i/v  and  x  =  uol  notation,  the  equation  for  the  libration 
points  (60)  becomes : 


3.2.2.  Outdance 

It  was  shown  in  section  3.1.1  that  the  number  of  guidance  functions  (ft 
equals  the  number  of  degrees  of  freedom.  Since  the  planar  restricted  three 
body  problem  has  two  degrees  of  freedom,  tpi  and  7)3  are  to  be  introduced 
and  evaluated.  Eqs.  (46)  in  their  modihed  forms  are: 


f  =  r(d  -f  n)2  -  +  f\(r,  B) 

1  hF 

rB=-  2r(6  +  ~ 


(63) 


Reduction  of  this  problem  to  a  second  order  differential  equation  is  possible 
only  if  the  Jacobi  integral  is  applicable,  that  is  if  and  973  are  components 
of  a  conservative  field,  or  if 

dfi  dipzr 

W  ^  “dT' 


In  other  words  if  the  guidance  functions  can  b*3  derived  from  a  potential, 
0,  then  the  potential  function  (F)  of  the  original  problem  can  be  modified 
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and  the  equations  of  motion  become 

6F* 

/ 

[in] 

(64) 

r0=  -2r{6^.  n)--^  \ 
r  oO 

where 

F*=F  +  0(r,  0). 

(65) 

This  way,  the  fourth  order  system  can  be  reduced  to  a  second  order 
differential  equation,  containing  r,  r',  r',  0,  0,  0f,  0,.  It  is  noted  that  by 
introducing  0  instead  of  using  inde|)eiKlent  qpi  and  functions,  the  problem 
is  restricted  to  conservative  gni<lance  forces,  but  this  restriction  is  not 
necessary  and  it  is  made  for  convenience’s  sake  only.  The  radial  and  tangential 
components  of  the  guidance  force  are  determined,  once  the  guidance  potential 
(0)  is  found,  by 

b0 

1  60 


The  .second  order  differential  equation  of  the  orbit  can  be  established 
simply  by  sub.stituting 

(9*)2  =  2(C-V*)^  1'  (07) 

and 

F*-  F  +  0 


in  eq.  (61).  Since  V(u,0)  is  given  (44)  and 


6u 


do* 

60 


;-^(Fe^0e), 


the  basic  guidance  equation  from  (61)  becomes 


g*{\i  )  M')j 


(68) 
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This  is  the  equation  which  corresponds  to  eq.  (18)  or  (19).  The  differential 
eq.  (68)  is  of  the  second  order  since  the  modifying  function  was  derived 
from  a  potential.  For  given  u  =  u(0)  this  equation  is  a  partial  differential 
equation  for  the  guidance  potential  (0).  The  erroneous  observation  might 
bo  made  that  partial  derivatives  would  not  appear  (seeeq.  (19))  if  a  potential 
function  for  the  guidance  forces  would  not  have  been  introduced.  On  the 
other  hand  it  should  be  observed  that  eq.  (18)  which  did  not  use  a  potential 
function  for  the  modifying  functions,  does  contain  partial  derivatives  and 
so  does  eq.  (25).  Therefore,  we  conclude  that  the  guidance  equation  is  always 
a  partial  differential  equation  either  for  the  guidance  functions  or  for  the 
guidance  potential. 

At  this  point,  to  illustrate  the  method,  we  select  an  orbit  and  show  how 
the  guidance  potential  is  determined.  Consider  the 

tt  =  *4-/9  cos  0  (69) 

conic  sections  and  require  that  the  mass  m  describe  this  orbit. 

It  can  be  shown  that  conic  sections  do  not  satisfy  the  differential  equations 
of  the  restricted  three  body  problem  [6],  therefore,  selecting  these  orbits, 
triviality  (0  =  0)  is  avoided.  The  practical  significance  of  such  orbits  is  also 
demonstrated  in  [6]. 

Substituting  eq.  (69)  into  (68)  the  partial  differential  equation  for  0 
becomes : 

60  60 

(70) 

From  this  equation  0(m,  0)  is  determined  keeping  in  mind  that  the  in¬ 
vestigation  refers  to  a  specific  type  of  orbits,  i.e.  to  those  which  are  described 
by  eq.  (69).  By  selecting  any  other  family  of  curves  for  the  path  and  sub¬ 
stituting  these  into  eq.  (68).  the  resulting  equation  will  have  the  same  general 
form  as  eq.  (70)  and  only  the  L(<P,  9)  and  M{0)  functions  will  change.  It  is 
to  be  noted  that  the  dependence  of  L  on  0  is  always  of  the  same  form,  i.e. 

L(0,  0)  =  P{9)  Q(9)  -  20  ^  <PR{9)  +  S(0)  (71) 

and  for  different  paths  only  the  P,  Q,  R,  S  functions  will  be  different. 

The  partial  differential  eq.  (70)  is  solved  by  the  standard  method,  i.e. 


dO 

^  d0 

(72) 

M{0) 

M<P.  0) 

'  ■''' 

M(0)  * 

(73) 

from  which 
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0  -i?(0)=!^2 

where  0  =  i?(O)  is  the  solution  of 

d0  L{0, 0) 
dO  “  M(0)  ' 

The  general  solution  of  the  partial  differential  eq.  (70)  is 


[m] 

(73) 

(74) 


or 


0  =  f2(O)  + 


K-/ 


_d0_\ 
i/(0)/  ’ 


(76) 


where  ^  is  an  arbitrary  function  of  its  argument. 

The  determination  of  the  actual  guidance  force  components  is  made  by 
means  of  eqs.  (66): 

,60 


and 


=  —  tt*  —  = 

Ou 


(76) 


For  the  conic  section  orbits  the  if,  P,  Q,  R,  S  functions  are  as  follows: 

^  sin  0 


if(0)  = 


P(0)=-2n 


+  cos  0)* 

1^  a*  -4-  +  2<xfi  cos  6 


and 


(«  +  P  cos  0)* 

Q{0)  =  g^0)  as  shown  in  eq.  (62), 

“  a*  +  /?»  +  2a/3  cos0 

S(0)  =  -(M  +  \QR). 

The  characteristic  eqs.  (73)  furnish  the  solutions: 

'Pi  =  u  —  ^  cos  0  —  2a  In  sin  0 - 3-^  In  tan  - 

P  ^ 


(77) 

(78) 


(79) 

(80) 

(81) 


and 


<^2  =  0-0(0) 
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where  0(0)  is  obtained  from  an  equation  corresj)onding  to  eq.  (74): 
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(10 

(10 


0_  +  _. 


(82) 


This  equation  can  l)e  reduced  to  an  Abel  ty|)e  differential  equation  and 
this  w’ay  the  guidance  pot  mtial  according  to  eq.  (75)  is  determinori,  subject 
to  the  selection  of  an  arb  trary  function. 


3.2.3.  Encke  method 

Three  body  trajectories  which  connect  the  ncighl)orhood  of  wii  with  the 
neighborhood  of  m2  are  of  considerable  practical  interest.  Such  orbits  often 
are  constnicted  of  three  parts.  The  first  j)art  is  a  two  body  aj)proximation 
fitted  to  that  part  of  the  trajectory  which  is  in  the  vicinity  of  m\.  The  third 
part  is  another  tw’o  body  a()proximation  fitted  in  the  neighborhood  of  m2. 
The  .second,  in-between  part  of  the  trajectory,  the  matching  of  those  two 
body  fits  and  the  j)roblem  of  approximating  the  trajectory  in  the  region 
where  the  force  field  is  e.s.sentially  a  three  body  field,  is  seldom  treated.  Since 
the  classical  Encke  methotl  or  the  variation  of  parameters  method  are  both 
ideally  suited  for  treating  domains  where  two  bo<ly  fiehls  dominate,  in  this 
chapter,  we  will  concentrate  on  the  “os.sentially”  throe  body  regime. 

Trajectories  which  in  a  rotating  coordinate  system  show  up  as  figure  8  or 
S-shaped  curves,  contain  intloction  points  and  in  these  regions  straight  line 
approximations  are  of  interest.  Neither  the  chussical  variation  of  j)arameters 
nor  the  conventional  Encke  method  is  designed  to  entertain  tho.so  regions. 

As  a  reference  orbit  the 


i’  =  a  sin  0  ^  ^  cos  0  (83) 

straight  line  is  used,  with  />  1 ll/n], 'X<(>  corros})onding  to  the  actual 
problem,  i.e.  the  line  will  intersect  the  ar-axis  (fig.  I)  between  the  origin  and 
m2  with  a  small  angle  of  jmsitivo  inclination.  Since  the  variable  u  is  re.serv'ed 
for  the  actual  orbit,  v  is  u.sed  for  the  reference  orbit  (comi)are  with  cqs.  (32) 
and  (33)).  Just  as  the  Encke  method  maintains  the  time  as  the  common 
indei>endent  variable  for  both  the  actual  and  the  reference  orbit,  in  the 
present  technique  the  indeiHJiident  angular  variable  will  tjo  coasidered  common 
to  both  orbits. 

Substituting  the  reference  orbit  into  eq.  (68),  one  finds  that  (since  t’-f  i’'  =  0) 
the  governing  equation  for  tl’.o  reforcnco  orbit  becomes: 


2n  \'(x^  1  =  (*  sin  0  -f-/3cos  0)^gt,*  —  [a/? cos  20  +  \('x'^  —  sin  20](Jq*.  (84) 
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The  general  solution  is  easily  obtained  by  the  method  of  characteristics ; 
g*  =  nv  +  P)  - 


where  is  an  arbitrary  function  and 


P  =  |/'a2+/?2  In  : 


lV  +  /?2-v 


2(p8in0‘- ACOflO) 

The  general  form  of  the  <p\  and  (p^  functions  require  now  the  determination 
of  0,  and  0,.  According  to  eqs.  (76): 

-vi0,  =  viig*g,*  +  Fr)  +  ^  / 

<Pt - v<Pt==v(g*g,*  +  Vt).  1 

This  completes  the  problem,  since  g*  is  given  by  (86)  and  V  by  (44). 
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